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The DC-response, namely the I-V and G-V charateristics, 
of a variety of composite materials are in general found to be 
nonlinear. We attempt to understand the generic nature of 
the response charactersistics and study the peculiarities as- 
sociated with them. Our approach is based on a simple and 
minimal model bond percolative network. We do simulate 
the resistor network with appropritate linear and nonlinear 
bonds and obtain macroscopic nonlinear response character- 
istics. We discuss the associated physics. An effective medium 
approximation (EMA) of the corresponding resistor network 
is also given. 

PACS numbers: 64.60. Fr, 64.60.-i, 05.70.Fh 



I. INTRODUCTION 

For measuring the susceptibility of a physical system, 
one generally studies its linear response for a small ex- 
ternal perturbation such that this perturbation does not 
appreciably change the basic nature of the system gov- 
erning the response characteristics under study. For an 
appropriately large perturbation (and a nondestructive 
one) the physical system may acquire newer modes of 
response and hence the system's response characteristic 
may change with the externally applied field. 

Composite systems, which are the object of our study 
here, may be thought to be comprised of many micro- 
scopic elements or grains (much larger than atomic di- 
mensions) having different physical properties. In the 
electrical case, it means that some of them may be metal- 
lic, some insulating and yet some others semiconducting. 
Further, the interface between them may also have differ- 
ent characteristics. If one takes the charge carrier (elec- 
tron or hole) as a random walker, a metal corresponds to 
diffusive motion giving rise to Ohm's law. This diffusive 
behavior changes in the presence of tunneling across a 
barrier, like M - 1 - M , S - M - S , M - s - M, s - / - s, 
etc. junctions {M stands for a metal, S a superconductor, 
and s a semiconductor). For example, a potential barrier 
called Schottky barrier arises within a semiconductor in 
intimate contact with a metal and gives rise to current 
i oc exp(u/i?o), where Eq is function of temperature. In 
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the case of a p-n junction heavily doped with Sb donors, 
the current seems to go from one linear region to an- 
other in steps, or the conductance smoothly connects two 
plateaus. Qualitatively similar curves |l| are observed 
for InP tunnel diode, Au doped Ge tunnel diode, etc. 
Further there could be inelastic processes by which an 
electron loses eneregy to the vibrational modes of species 
within the barrier. These processes thus provide extra 
paths or channels for current flow (thus increasing the 
conductance), and hence one may observe breaks in i-v 
curves joining two piece-wise linear (ohmic) regions. In 
the following percolative approach to obtain the direct 
current (DC) response of composite systems, one treats 
these elements semi-classically in the sense that one does 
not solve for the wave-function but takes the quantum 
mechanical effects indirectly by including only the tun- 
neling currents through barriers 

It is clear that if these microscopic tunneling ('nonlin- 
ear') elements become linear as discussed above beyond 
some characteristic microscopic voltage, then the nonlin- 
ear conductance of the whole system made up of such ele- 
ments would change over to ohmic behavior (constant dif- 
ferential conductance) beyond some system-specific ap- 
plied voltage. The point which may not be apparently 
obvious is that even if the nonlinear bonds remain so all 
the way upto infinity (i.e., their resistances become zero 
('super'-conductor) at infinite voltage or frequency), the 
macroscopic system comprising of such bonds still satu- 
rates at high voltage to another linear region or ohmic 
behavior if the tunneling bonds do not percolate by them- 
selves. In either case, the resistive behavior of the origi- 
nal percolative structure ( zero voltage ) guides the resistive 
behavior of the 'shorted' percolative structure. Thus, for 
vanishing driving fields, the system is in some state pro- 
viding a fixed number of channels for the response, but 
asymptotically for infinitely large driving fields, the sys- 
tem is in another state representing a much higher num- 
ber of channels. The generalized susceptibility (which 
could be thermal or electrical conductivity, any elastic 
moduli, or permeability, etc.) of this two (or, multi)-level 
system crosses over from one asymptotic linear region 
(value) to another (of higher saturation value) through 
some nonlinear region. It may be noted here that in 
rupture type of breakdown {e.g., electrical fuse 1^) phe- 
nomena, essentially the opposite thing happens, namely 
the system crosses over from a higher susceptibility value 
(highly efficient and stable network) to a lower one where 
the system cannot hold itself anymore. In electrical ter- 
minology, the system as a whole fuses to give rise to an 
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insulator, or in mechanical terminology, the structure be- 
comes mechanically unstable. Thus this enhancement 
or (de- enhancement) of the number of system- spanning 
channels and their eventual saturation as a function of 
some driving field is at the heart of the nonlinear response 
which is considered here in this work. 



A. Experimental Facts 

Nonlinear transport characteristics have been revealed 
in many different types of materials, including an early 
work by van Beek and van Pul (1964) on carbon-black 
loaded rubbers. References in this regard are meant 
to be representative and not exhaustive. We discuss 
below some interesting and common features of com- 
posite materials, specially which are highly structured 
and give rise to some sort of universal behavior. Com- 
posite materials are typically a mixture of two or more 
phases, namely a mixture of metallic kind of material 
in an insulator where the mixture is not in the atomic 
scale. The metallic islands formed in the insulating 
matrix (as may be found by tunneling electron micro- 
graph, TEM) are typically of dimensions much greater 
than the atomic size but much smaller than the macro- 
scopic system size. As an example we may here re- 
fer to the carbon-black-polyvinylchloride (PVC) compos- 
ites 1^. Carbon blacks composed of small but compli- 
cated shaped particles usually exist in the form of "high- 
structure" aggregates, whereas smaller and geometrically 
simpler particle aggregates are also possible and they are 
called "low-structure" blacks. A schematic description 
of such composites made of different forms or structures 
are given in Ref. [^). The conductivity exponents for 
those three types of carbon-black-PVC samples were also 
measured. They were found to be t = 4, i = 2.8 and 
t — 2 for "low-structure" (commercially called Mogul- 
L black), "intermediate-structure" (Cabot black) and 
"high-structure" composites respectively (8|. Only the 
"high-structure" composite [t = 2) was found to be in 
the universality class of ordinary percolation problems 
and we actually have those kind of systems in our mind. 
In fact the following features pertain to a wide variety of 
such composite systems available. 

• Very Low percolation threshold: Usually these 
composites exhibit an unusually low percolation 
threshold. For example, in an experiment on 
carbon- wax system Q , there is a percolation tran- 
sition at Pc = 0.0076. Very low thresholds are 
also reported for other systems, e.g., Pc — 0.002 
for carbon-black-polymer composites |^ and Pc = 
0.003 for sulphonated (doped) polyaniline networks 

• Qualitatively identical I-V (as well as dl/dV 
vs. V) response both below and above 



threshold: As an example one may note the ex- 
periment on Ag particles in KCl matrix by Chen 
and Johnson [Q, where similar nonlinear transport 
was reported both below and above pc ~ 0.213, 
even though the exponents for these two regions 
were reported to be grossly different. The fact 
that charge is carried even below pc indicates that 
tunneling/ hopping between disconnected conduct- 
ing regions does take place and is responsible for 
the nonlinearity. 

• Power-law growth of conductance: For small 
V, the excess conductance (above the constant/ 
ohmic part) is claimed to grow as a non-integer 
power law with exponent 6 = 1.36 in 3D 
whereas previous theoretical works took 5 to be an 
integer equal to 2 to keep voltage inversion sym- 
metry. Power-law in conductance against voltage 
{G-V) automatically implies the power-law in I- 
V characteristics. The nonlinearity in I-V curves 
and the associated power-law have been observed 
in a variety of experiments (although in most of 
the experiments the G-V curves are not examined). 
Rimberg et ai, ||ll| measured the I-V characteris- 
tics of one- and two-dimensional arrays of normal 
metal islands connected by small tunnel junctions 
and found a power-law: I ^ {V — Vg)", Vg being 
the system's threshold voltage. The nonlinearity 
exponent a is found to be 1.36 ± 0.16 in the case 
of ID array and 1.80 ± 0.16 for 2D. However, the 
values of a were earlier predicted theoretically to 
be 1 and 5/3 in ID and 2D respectively by Middle- 
ton et al. [ p^ , through an approximate analytical 
calculation. 

• Saturation of conductance: In general for com- 
posite systems and for many other materials the 
G-V characteristics show a very interesting behav- 
ior. The G-V curve is seen to saturate for an ap- 
propriately high enough voltage below the Joule- 
heating regime. The typical curve then looks like 
a nonlinear sigmoidal type function interpolating 
two linear regimes. Some of the experiments 

as mentioned above present the G-V data in this 
form for the purpose of proper analysis. A similar 
nonlinear curve for conductivity against field is also 
presented by Aertsens et al. [ [l3| for microemulsions 
of water in oil under an external electric field. 

• Crossover current for nonlinearity: The cur- 
rent Ic (voltage Vc) at which the conductance G 
differs from the zero-current (zero- voltage) conduc- 
tance Go by some small 1%) but arbitrarily cho- 
sen fraction e is called the crossover current (volt- 
age). It is seen to scale as Ic oc Gg (K oc Gq~^), 
where x is called the crossover exponent for non- 
linearity. For carbon- wax composites, x = 1.4 

in 3D, and for discontinuous thick gold films, x = 
1.5 as found by Gefen et al. |4|. One can easily 
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check that this exponent is related to S above by 
<5 = l/(x- 1). 

• Temperature-dependent conduc- 
tion: Composite systems display very interesting 
temperature-dependent conduction properties par- 
ticularly in the low temperature regime where the 
conduction is mainly due to phonon-assisted hop- 
ping of the electrons between randomly spaced lo- 
calized states. One then needs to consider Mott's 
variable range hopping (VRH) or some of its varia- 
tions namely log G ~ T~'^ for low T. The conduc- 
tance of the sample goes through a maximum as 
the temperature is increased still further towards 
some type of metallic behavior. 

The first thing to note from the enlisted facts is the 
ultra-low percolation threshold even when the included 
conducting phases are isotropic. This alongwith the fact 
that many of these nonlinear systems carry current even 
below pc indicates strongly that tunneling through dis- 
connected (dispersed) metallic regions must give some 
virtually connected percolating clusters. ^From the non- 
linear I-V characteristics (e.g., see the experiment by 
Chen and Johnson Q) it is observed that the response 
behavior is reversible with respect to the applied field 
in the sense that the response curve (current or voltage) 
does almost trace back as the field is decreased. This fact 
also indicates that reversible tunneling is responsible for 
such a behavior. Also the temperature-dependent con- 
ductance with a maximum at some characteristic temper- 
ature (dependent on the amount of disorder present) and 
the Mott VRH type behavior at very low temperatures 
give further credence to tunneling assisted percolation. 
So, in the following section, we propose a semi-classical 
(or, semi-quantum) model of percolation which works 
on the borderline between a classical and a quantum pic- 
ture. Quantum physics enters our discussion through the 
possibility of tunneling of a charge carrier through a bar- 
rier (which docs not exist classically). 

One may consider a wide variety of nonlinear transport 
(apart from the nonlinear electrical conduction) processes 
and study them under the framework of percolation the- 
ory. One also exploits the key features and ideas of one 
such model and incorporate them into the other. As 
an example, one draws analogy between laminar flow in 
tubes and the electrical currents and uses the language of 
electrical network for convenience. The volumetric flow 
rate q is identified with the current i and the pressure 
drop AP across a tube or a pore to the voltage difference 
V across a bond. Thus, for example, the flow of polymers 
is modelled (see e.g., Ref. and references therein) 
considering a power law i — gv°' to each bond (tube or 
pore), where g is a generalized conductance. Consider- 
ing the percolating network model of porous media one 
employs Monte-Carlo simulations or an effective medium 
approximation to calculate the rheological properties of 
a power-law fluid in flow through porous media. Foam 



is a non-Newtonian fluid which is used in displacement 
and enhancement of oil recovery from the porous rocks. 
However, to move the foams through the pores, external 
pressure has to exceed a certain critical value P^. In brit- 
tle fracture, no microcrack nucleation process takes place 
unless the applied stress exceeds a critical value for the 
system. A convenient description of this particular prob- 
lem is done through the electrical analogoue of breaking 
with the introduction of random fuse model |^ . One de- 
fines the 'fuse' as a device with a constant conductance 
when the applied voltage across it is less than a certain 
critical value Uc, beyond which it is an insulator. This 
electrical network model is a scalar analogue of the cor- 
responding vector elasticity problem where the former is 
simpler to handle. 

In most of the works we mentioned above, dilution 
plays an important role. These are models where the 
percolation theory has been the underlying framework 
and many of the interesting properties may be related 
to the cluster statistics. In this paper we first discuss 
briefly our proposed model for studying the effective non- 
linear conduction properties of various composites. We 
then present the nonlinear current-voltage {I-V) charac- 
teristics in our model system and discuss the associated 
details. 



II. THE MODEL AND ITS PERCOLATIVE 
ASPECTS 

To mimic charge transport in composite systems or dis- 
persed metals, we assume the grains (metallic or metal- 
like) to be much larger than atoms but much smaller 
than laboratory-scale macroscopic objects. These grains 
are randomly placed in the host material which are in- 
sulators. Further we assume quantum mechanical tun- 
neling between such grains across some potential bar- 
riers. Clearly these barriers would depend on the lo- 
cal geometry of the insulating and the metallic grains. 
Since in practice, the tunneling conductance should fall 
off exponentially, the tunneling should have some length 
scale designating an upper cut-off (for tunneling to oc- 
cur) in the separation between two metallic grains. Fur- 
ther, the separation between two grains or the poten- 
tial barrier can vary continuously between zero and some 
upper cut-off. For simplicity and to capture the basic 
physics, we construct a bond (lattice) percolation model 
for this problem, such that tunneling may take place only 
between two nearest neighbor ohmic conductors (or, o- 
bonds) and no further. Thus one may imagine a virtual 
bond sitting at each such gap which conducts current 
nonlinearly due to tunneling phenomenon. We call these 
tunneling conductors as tunneling bonds (or, t-bonds). 

We stress again that our approach would be to solve an 
appropriate electrical network based on a semi-classical 
(or, semi-quantum) percolation model. Made of both 
random resistive and tunneling elements, this network 
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would be called a random resistor cum tunneling-bond 
network (RRTN) Now tunneling may take place 

through the tunneling bonds in various ways, so that 
the functional form of the tunneling current as a non- 
linear function of the potential difference across them 
may be quite complicated. For simplicity, we address the 
aspects of nonlinearity in a macroscopic system which 
comes through two piecewise linear regions of a t-bond 
such that the conductance is zero upto a threshold volt- 
age and a given non-zero value beyond. In this context, 
it may be noted that disorder in many quantum systems 
like charge-density-wave (CDW) systems or flux-vortex 
lattices of type-II superconductors can give rise to 'pin- 
ning' or inhibition to transport upto a critical value of 
the applied field above which tunneling is active. The 
piecewise linear transport is in fact a highly nonlinear 
process as there is a cusp singularity at the intersec- 
tion point. The transport due to tunneling which is the 
source of nonlinearity in the experimental systems [^|J^] 
we focus on, can be well approximated in this way and 
thus the nonlinearity of the macroscopic systems may be 
understood at a qualitative (and, sometimes even at a 
quantitative) level. Next, one notes as discussed above 
that in many physical systems, the response is negligibly 
low (or there is no response at all) until and unless the 
driving force exceeds a certain threshold value. A class 
of problems exist where sharp thresholds to transport 
occur. The examples in the electrical case are a Zener 
diode, a CDW system or a type-II superconductor and 
in the fluid permeability problems, for example, a Bing- 
ham fluid (where there is a critical shear stress Tc, above 
which it has a finite viscosity and below which it is so 
enormously viscous that it does not flow). In our RRTN 
model, we work with i-bonds which have zero conduc- 
tance below a threshold; see Fig. |l|. 



FIG. 1. The piecewise linear current-voltage (i-v) charac- 
teristic of each tunneling conductor. 

Let us look at the uncorrelated random bond percola- 
tion first. In this model, the conducting o-bonds are oc- 
cupied randomly with a certain volume fraction p. One 
uses a random number generator which generates ran- 
dom numbers distributed uniformly in the interval [0, 1], 
to fix the positions of o-bonds in a certain configuration. 
So the system now has p fraction of occupied bonds and q 
(= 1-p) fraction of unoccupied bonds. Upto this it is the 
standard (uncorrelated) bond percolation problem where 



the percolation threshold in a square lattice is Pc = 0.5. 
A geometrical connection (or a system-spanning cluster 
of o-bonds) is established from one end of the 'config- 
uration averaged' system to the other only at or above 
this volume fraction in an infinite size system. If the oc- 
cupied bonds are ohmic conductors and the unoccupied 
bonds correspond to insulators (the host) then we have 
the corresponding random resistor network (RRN) which 
carries an average non-zero current only when the volume 
fraction, p > pc- Clearly this response is linear or ohmic. 

But, our percolative model is not just a random mix- 
ture of two phases. For our convenience we take a square 
lattice in 2D. The basic physics should remain the same 
if we go over to 3D. As usual, the o-bonds are thrown 
at random at a certain volume fraction p. The rest 
{\ — p) fraction contains insulators. Now we allow tun- 
neling bonds (t-bonds) only across the nearest-neighbor 
{nn) gaps of two conducting bonds (and no further) if 
an appropriate voltage is applied externally across two 
opposite sides of the RRTN. The RRTN model is now 
comprised of ohmic conductors, pure insulators and some 
non-ohmic (tunneling) conductors. To have a clearer 
view of the proposed model system we refer to the Figs. ^ 
where we show all the ohmic bonds thrown at a concen- 
tration p, and all possible t-bonds. In the Fig. ||(a) , we 
show a realization of a square lattice at p = 0.15 in 2D 
where the system does not have any connectivity from 
one end to the other even if we consider all the i-bonds 
to be active. In the Fig. ||(6) aX p = 0.25, there is a sys- 
tem spanning cluster (network) of o-bonds and i-bonds 
even though the system is actually below the geometrical 
threshold {p < pc). This demonstrates that the effective 
percolation threshold of the system is lowered as the 
t-bonds come into play. Thus even if the system is be- 
low Pc, and hence has zero conductivity at a vanishingly 
small voltage, its conductivity will be non-zero and keep 
growing at appropriately higher voltages, and hence the 
system will behave non-ohmically. If the system is al- 
ready above the geometrical threshold (pc = 0.5) then 
the i-bonds simply add extra paths to the already exist- 
ing current carrying network (backbone) of o-bonds and 
give rise to similar non-ohmic behavior. 

If we consider the percolative aspects of the model sys- 
tem at very large voltages (suppose all the t-bonds over- 
come their threshold and conduct ohmically), then this 
is no more a pure bond percolation problem of a random 
binary mixture. Rather we may think of it as a very spe- 
cific correlated bond percolation problem. This is because 
the positions of the t-bonds are totally correlated to the 
positions of randomly thrown o-bonds. The disorder is 
in the position of the o-bonds only. Once the positions of 
the o-bonds are given for a particular configuration, the 
positions of the t-bonds are automatically determined for 
that configuration in our model. We have found out ear- 
lier 1 18 1 that the new percolation threshold with all the 
<-bonds added is Pct ~ 0.181. We have also addressed |T^ ] 
the question of universality class of this problem in this 
limit and found that the problem belongs to the same 
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universality class as that of the pure (uncorrelated) bond 
or site percolation. 




FIG. 2. Tow typical configurations of a square lattice of 
size 10 X 10 with p < pc and an arbitrary but very large 
voltage where all the tunneling bonds indicated by broken 
Unes are active, (a) For p = 0.15 the system is seen to 
have no connecting path across its two ends even when one 
considers the tunneling bonds. (6) This configuration at p 
= 0.25 is seen to have no connecting path with only ohmic 
bonds, but has at least one such path comprised of both the 
ohmic and the tunneling bonds. 



III. EFFECTIVE MEDIUM APPROXIMATION 
(EMA) FOR THE CORRELATED BOND 
PERCOLATION MODEL 

An effective medium approximation gives reasonable 
result for any p away from the threshold (here Pct)- The 
EMA has been used to calculate fairly accurate conduc- 
tivity behavior for general binary mixtures except in the 
vicinity of the critical regions. This approach is old and 
had been devised for the transport properties of inho- 
mogeneous materials first by Bruggcman p9| ] and then 
independently by Landauer [ pO[ . Its successful applica- 
tion to the percolation theory by Kirkpatrick has 
drawn attention of many others in this field. It has then 
been applied for a wide variety of inhomogeneous materi- 
als. The development has some similarity with the well- 
known coherent potential approximation (CPA) for treat- 
ing the electronic properties of the binary alloy problem. 
The composite systems have seen a series of attempts in 
this direction while dealing with the nonlinearity in the 
response involved (see Ref. |^ and references therein). 

Here we follow the method as described in ref. pl[ . 
The method is as follows. Consider a random electri- 
cal network on a hypercubic lattice of dimension d > \ 
{d ~ 2 for 2D, d = 3 for 3D etc.). The basic idea is 
to replace a random network by a homogeneous effective 
network or an effective medium where each bond has the 
same average or effective conductivity Gg. The value of 
the unknown Ge is calculated in a self-consistent man- 
ner. To accomplish this, one bond embedded in the ef- 
fective medium is assigned the conductivity distribution 
of the actual random network. The value of Ge is then 
determined with the condition that the voltage fluctua- 
tion across the special bond within the effective medium. 



when averaged over the proper conductivity distribution, 
is zero. The voltage (v) developed across such a special 
bond can be calculated |^ for a discrete lattice of z (= 
2d for a hypercubic lattice) nearest neighbors as 



vx{Ge-g)/[g + (d-l)Ge 



(3.1) 



The requirement is that the average of v be zero when 
the conductance for the special bond may take any of the 
ohmic, the tunneling or the insulating bond value with 
appropriate probabilities for the actual network. 

The probability of a bond to be ohmic, tunneling or 
purely insulating according to the considerations of our 
model is given below for 2D square lattice. 



Po=P, 

Pt = (p^ + 3p\ + 3pq^fq, 

P, = [l-{p'' + 3p\ + 3pqY]q, 



(3.2) 
(3.3) 
(3.4) 



where q = 1 — p. For a distribution of three types of 
resistors, we have 

fig) = PoS{g - go) + Pt5{g - gt) + P^5{g - g^), (3.5) 

where (^o, gt and gi are the conductances of ohmic, tunnel- 
ing and the insulating resistors respectively. The EMA 
condition stated above i.e., < v > — Q reads 



dgf{g){Ge^g)/[g+{d~l)Ge]^Q. 



(3.6) 



Putting Eq. ( pT5| ) into the above Eq. (3.6), we get 

Po{Ge~go) , PtiGe-gt) , P^{Ge-g^ 



[go + {d- l)Ge] [gt + {d - l)Ge] [g^ + {d - l)Ge] 

(3.7) 

The above equation reduces to the quadratic equation 

AGl + BGe + C ^0, (3.8) 

where A = {d - 1)^, B ^ [{d - l){Pt + Pi)go + {d - 
l){Po + P,)gt -{d- l)\Pogo + Ptgt)] and C = [Pm - 
{d — l)gogt{Po + Pt)], considering the conductance of the 
insul atin g bonds to he gi — 0. The valid solution of 
Eq. (ph is 



G„ 



B + {B^ - 4AC) 
2A 



1/2 



(3.9) 



Now one can obtain the linear conductance of the macro- 
scopic model composite system in 2D and 3D, putting 
d=2 and 3 respectively, given some specific values or 
functional forms for the conductances of the elementary 
components like the ohmic and the tunneling bonds. 

In our model we have assumed that the conductance of 
a tunneling bond, when it overcomes its voltage thresh- 
old, is the same as that of an ohmic bond. We believe 
that this assumption does not change the phase transi- 
tion characteristics. Now in the limit of large enough 
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external voltage if we assume gt = .go = 1 and gi = 0, we 
have 



d{Po + Pt)-l 
(d-l) ■ 



(3.10) 



At the percolation threshold pd the effective conductance 
(Ge) of the system is zero. In 2D we now have the fol- 
lowing equation putting d = 2 in Eq. (3.10), 



2{Po + Pt) = 1. 



(3.11) 



Solution of the above equation gives p = Pct =1/4. This 
may be compared with our simulation result pq ] which 
gives Pct — 0.181. The EM A is basically a mean field 
calculation which overestimates the percolation thresh- 
old value in lower dimensions. However, for pure bond 
percolation in 2D, this gives the exact result Pc = 1/2. 
This is so because the square lattice in 2D is self-dual in 
the case of purely random bond percolation problem. 

We may also find the value of Pct for 3D where the 
probability of tunneling bond is 



Pt = ip' 



5p*q 



10/ (7^ 



lOp'^q^ 



^pq^q 



(3.12) 



Usin g d = 3 in Eq. (3.10), an equation similar to 
Eq. ( 3.1l| ) is obtained, which when solved gives pct = 1/8 
exactly. 



IV. THE NONLINEAR I-V CHARACTERISTICS 

For the work presented below, we make the simplify- 
ing assumption that all the tunneling bonds (t-bonds) in 
our RRTN model have an identical voltage threshold (vg) 
below which they are perfect insulators and above which 
they behave as the ohmic bonds. One could certainly 
introduce disorder by making Vg random as in Ref. [ p3| . 
In our case disorder is already introduced through ran- 
dom positioning of the bonds, and we believe that our 
assumption should not affect the dilution-induced non- 
linearity exponents. Indeed, as shown in the sequel, our 
model gives richer possibilities (dilution dependence) for 
the nonlinearity exponent in a composite system. 

As the externally applied field (DC) is increased be- 
yond some macroscopic threshold, some of the i-bonds 
overcome their microscopic thresholds and may thus in- 
crease the overall conductance of the system if the pro- 
cess leads to newer parallel connectivities for the whole 
macroscopic composite. Our computer simulation to 
study these effects |24 involves the solution of Kirch- 
hoff's law of current conservation at the nodes of the 
RRTN with the linear and the nonlinear (assumed piece- 
wise linear) resistors and the standard Gauss-Seidel re- 
laxation. Current response was averaged over 50 con- 
figurations in all the cases mentioned here. We obtain 
current (/) against voltage V and therefrom the differen- 
tial conductance (G = dl/dV) for the whole network at 



a given volume fraction p of the ohmic bonds. Simula- 
tion results for nonlinear I-V curves for a square lattice 
of size L — AO are plotted in Fig. ^ for p — 0.3 to 0.9. 
For p > Pc, one may note that the I-V curve is linear up 
to a certain voltage (Vg), beyond which the nonlinearity 
shows up. For p < pc, there is no current (zero con- 
ductance) below a threshold voltage (Vg), beyond which 
the nonlinear conduction starts. Nonlinearity is always 
there in the I-V response for any value of p in the in- 
terval Pct < p < 1. However, for pct < p < Pc, there 
is no system-spanning path with the ohmic bonds (in an 
average sense) and the response (average current) is zero 
for small voltages. The response in this case starts out 
nonlinearly from a non-zero average threshold voltage. 
On the other hand, for p > pc, the system always has a 
conducting path through the ohmic bonds (again on an 
average), and so up to a certain voltage (Vg) there is a 
constant non-zero conductance and the average response 
is linear. As V is increased beyond Vg, more and more 
current carrying paths are added with the help of the 
i-bonds and in effect the conductance starts increasing 
with V. 
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FIG. 3. A set of nonlinear I-V curves for p = 0.3 to 0.9 for 
the RRTN network of size 40 x 40. The current response is 
averaged over 50 configurations each. 

We find that the entire nonlinear regime of a I-V curve 
can not be fitted by a simple power-law which in general 
may be fitted by a polynomial function. Even after doing 
that, the exponent of nonlinearity from the fitting of the 
various I-V curves remained somewhat ambiguous since 
the fitting was not very robust. Thus to have a better 
idea, we obtain the differential conductance (G = dl/dV) 
by taking the numerical derivative of the the I-V data. 
In these derived G-V curves, one can identify the thresh- 
old voltage Vg as the onset of nonlinearity and the onset 
of the saturation regime (G/) much more clearly than in 
the I-V curves. The G-V curves corresponding to p — 
0.7 and 0.3 in the Figure ^ are shown in the Figs. ^ and 
H respectively. As expected for the two different cases, 
in Fig. m (p > Pc) the initial conductance has a finite 
nonzero value (Go ^ 0), while in the Fig. ^ (p < pc) 
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the initial conductance is zero (Go = 0). The general 
nature of the G-V characteristics is that eventually all 
of them become flat, i.e., assume a configuration- and p- 
dependent constant conductance beyond some very large 
voltage Vs. The reason is that for a finite-sized system 
there is always a very large but finite voltage Vs at which 
the conductance of the whole system saturates to an up- 
per maximum value G/ (another linear regime) since all 
the possible t-bonds become activated and no more chan- 
nel parallel to the backbone come into play for a further 
increment of V. Below we concentrate on the analysis of 
the nonlinear conductance behavior of the model system 
in an effort to first find a general functional form to de- 
scribe all these different G-V curves and then to find the 
exponent of nonlinearity at different volume fractions. 
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FIG. 4. The G-V curve corresponding to p = 0.7 of Fig. || 
(p > pc). The threshold voltage iVg) is indicated in the Fig- 
ure. Initial conductance (Go 7^ 0) and the final conductance 
(G/) are also indicated. 




FIG. 5. The G-V curve corresponding to p = 0.3 of Fig. |^ 
(p < Pc). The threshold voltage (Vg) is indicated in the Fig- 
ure. Initial conductance Go = and the final conductance 
(G/) is as indicated. 



A. Analysis of the Behaviour of Conductance 

Guided by the conductance behavior of a complex sys- 
tem made out of many simple prototype circuits jlj] and 
by the fact that the initial power-law law type growth of 
G beyond Go finally saturates to G/, we try to fit the 
whole region in our simulation by a function of the form: 



G — Go 

— Gf ~ Gd[^ 



iv > Vg), 



(4.1) 
(4.2) 



where AV ^ V — Vg is the shifted voltage from where 
nonlinearity appears. For concreteness, we discuss here 
the fitting of a sample data set for L = 40, p = 0.6. For 
this sample Go = 0.154, G/ = 0.881. The parameter for 
the best fit in this case as obtained by a simplex search 
procedure are A ^ 3.27 x 10"-*, fi = 1.408, and 7 = 125. 
Such large values of 7 are obtained for all the cases stud- 
ied and yet the approach to G/ was slower than in ac- 
tual data. This obviously indicates that the approach to 
saturation is not a power-law type and probably an expo- 
nential function is involved. At this point we refer to the 
Fig. ^ where we have shown a G-V curve for p = 0.8 and 
L = 40 for which the conductance (G) seems to saturate 
(to the naked eye) at a voltage above V = 20 and one 
can see a practically fiat regime in between y = 20 to 40. 
However, in the same figure we have demonstrated in the 
inset by zooming in on the y-axis that the actual satura- 
tion is yet to come and that the conductance (G) in this 
regime increases very slowly with V. The reason is that 
there are some tunneling bonds (typically in the trans- 
verse direction to the electric field) which do not become 
active even with the application of a very large voltage 
implying that the conductance for the whole system is 
yet to reach the complete saturation. This pathology 
supports the fact that the final saturation occurs at an 
extremely large voltage, and that the approach to the 
saturation is indeed very slow. 
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FIG. 6. A typical G-V curve for L = 40 and p = 0.8 with 
an apparent onset of the saturation regime beyond V « 20. 
The inset shows that the G does not yet saturate in the true 
sense and instead ever increases in the regime 20 < V < 40. 
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FIG. 7. The demonstr ation of the fitting of a G-l^ curve 
by the proposed function (4.7) shown for L = 40 and p = 0.8. 
The fitting is better with 7 7^ 1 (shown in dashed line) for the 
entire data set than that with 7=1 (shown in full line). 

We tried the following four plausible functions involv- 
ing exponentials for the entire nonlinear regime replacing 
the above Eq. (O). 



G 



Go + Gdexp(-A/AF) 

Go + Gdtanh(AAF") 

Gf - Gd exp[l - exp(AAl/")] 

Go + Gd[l-exp(-AAF'')]'' 



(4.3) 
(4.4) 
(4.5) 
(4.6) 



Out of all the opti mally fitted functions as described 
above, the last Eq. (4.6) is the best in the sense that 
it gives the minimum mean square deviations (MSD) for 
the G-V data. It may be noted that a special case of 



Eq. ( W ) with 7 = 1 was the form of nonlinearity used 
by Chen and Johnson in fitting their experimental 
conductance against voltage data for a composite system 
of Ag-KGl. But for all p and L's considered by us, 7 = 
1 was found inadequate for fitting the typical sigmoidal 
curves. For example we have shown in Fig. ^ the fit- 
ting with the function (46) for p = 0.8, L = 40: the 
restricted case of 7 = 1 shown by full line, and an unre- 
stricted, optimally fit 7 7^ 1 by dashed line. Clearly the 
unrestricted case fits the data extremely well and gives 
an MSD which is much smaller than that of the restricted 
case as mentioned above. 



B. The Nonlinearity and the Crossover Exponent 

The conductance (G) for metal-insulator composite 
systems starts growing nonlinearly with the applied volt- 
age V from an initial value (Gg) and finally saturates to 
a value G/ (as indicated in the Figures || and ^. Go is the 
conductance when none of the i-bonds is active; this hap- 
pens in the usual percolation model, without tunneling. 
Likewise, G/ is the conductance when all the so-called 
tunneling bonds (<-bonds) are actually active and taking 
part in the conduction. There are three distinct regimes 



which, in general, can be precisely located from the G-V 
characteristics than from the I-V charactersistics: 

• (i) Upto some voltage Vg the initial conductance 
(Go) of the system is either zero or a fixed finite 
value depending on whether the system has initially 
conducting path through the ohmic (or metallic) 
bonds or not (see Figs. ^ and ||). 

• (ii) Beyond Vg the nonlinearity starts showing up. 
The conductance (G) is nonlinear in the regime 

Vg<V< V,. 

• (iii) Beyond the voltage Vg , G saturates to a voltage 
G / (see Figs. |4| and |) . 

The data obtained through our model system in 2D were 
fitted through the following general formula as discussed 
above: 



G = Go + Gd[l - exp(-AAF'')] 



(4.7) 



where Gd ^ G f — Gq. Clearly, irrespective of the value 
of Vg, Go is the conductance in the limit ^ — > 0. Experi- 
mentally G / may be obtained by applying a large enough 
voltage such that Joule heating remains unimportant. In 
our computer simulation on finite sized systems, we find 
Vs to be a large but finite voltage (in arbitrary units). 
For example for squares with L=40, Vg is found to be 
typically of the order of 10''-10^. 

In a recent experiment by Chen and Johnson 0| in 
Ag-KGl composite (silver particles in KGl matrix), very 
similar G-V curves were obtained and the non-ohmic ef- 
fect was postulated to arise from a localized reversible 
dielectric breakdown between narrowly separated metal 
clusters in the metal-insulator composite. The interclus- 
ter or interparticle spacing may have some distribution 
which is related to the fractal dimension d/ of the net- 
work at the threshold. Chen and Johnson used the 
following conductance behavior for their data: 



G = Go + {Gf- Go)[l - cxp{-V/Vg) 



n{df) 



(4.8) 



which is a sp ecia l form of the function ( [4. 7] ) we propose. 
We find Eq. (^) to be inadequate for the representation 
of our and experimental data (see above). The exponent 
n{df) is in fact the same as S and is hence related to the 
nonlinearity exponent a (see Sec. |I A[ ) by a = n{df) + 1. 
It has been further shown in the above work that n{df) 
increases as the silver volume fraction is decreased and 
it shows a sharp change at the threshold. We shall dis- 
cuss in this section how our model captures this dilution- 
dependent nonlinearity exponent. 

For a meaningful comparison of all the G-V data with 
different Go, G/, Vg, etc., we look at the scaled con- 
ductance G — {G — Go)/Gd against the scaled voltage 
V = {V -Vg)/Vg. To see if the G-V data for various 
values of p both below and above Pc scale, we first looked 
within the range 0.48 < p < 0.52 {i.e., very close to pc), 
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and found that all the data do reasonably collapse. In 
the Fig. ||we show such a plot for a 20 x 20 system. This 
suggets the following general form for the functional be- 
havior; 

G = f{V), (4.9) 

where f{x) is a function such that /(O) = 0, and /(oo) = 
1 and is otherwise quite general as long as it represents 
the behavior of G very well. Clearly the scaled function 
in Eq. (4.7) satisfy these properties very well. 
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Thus the nonlinearity exponent 6 is related to the fitting 
parameters /i and 7 by 5 = nj. In our earlier work 
Jl^ with our preliminary observations we had reported 
the nonlinearity exponent S to be close to 1 and to be 
independent of p near the geometrical percolation point 
Pc = 0.5. Indeed, in most of the experiments an average 
value of the above exponent is reported for the data for 
samples close to Pc- 

Further careful analysis of the results at widely differ- 
ent volume fractions indicate that the nonlinearity expo- 
nent S increases significantly as we go sufficiently away 
from the percolation threshold (both below and above). 
This becomes apparent from the shape of the G-V curves 
for different volume fractions in a wide range of p (from 
0.3 to 0.9) in Fig. ^ corresponding to the I-V curves 
shown in Fig. ||. In Fig. |l^ we plot the scaled conduc- 
tance (G) against scaled voltage {V) corresponding to all 
the G-V curves in Fig. ^. The scaled data for all the 
curves now do not fall on top of each other indicating 
that all of them can not be described by the same fitting 
parameters /x and 7, even though the form of the best 
fitting function f{x) remains the same. Hence the non- 
linearity exponent in the power-law regime (AT^ 0"'") 
for these curves of different p are not identical. 



FIG. 8. The plot of scaled conductance against scaled volt- 
age for various p around pc- 
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FIG. 10. The scaled conductance against the scaled voltage 
for different p (both below and above pc) corresponding to the 
curves in Fig. pi The data do not coUapsse at all indicating 
different exponents for different p's. 



FIG. 9. The family of G-V curves corresponding to the 
Fig. |. 

Here we point out that the threshold voltage Vg is the 
only relevant variable that enters into the scaling func- 
tion. The other voltage scale Vs which is the onset voltage 
for s atur ation, is seen to have no role in the above scaling 



Eq. (|4.9| ). For small AV = V — Vg, i.e., near the onset 
of nonlinearity, the excess conductance AG = G — Go 
varies with the voltage difference {AV) as a power-law 
as one may easily check by expanding Eq. (4.7) in the 
limit (Ay 0): 

AV'^ = AV^. (4.10) 



AG 



The fitting of the individual G-V curves for squares of 
sizes L = 20 and 40 and some data for L = 60 and 80 
for different p in t he range of p = 0.2 to 0.9 were done 
by using Eq. (T7). The Fig. |ll| demonstrates that the 
exponent 6 for L = 40 increases from a value close to 1 
at p = Pc to values close to 3 and above on either side 
(0.3 < p < 0.9). We found that the value of 5 at pc lies 
in the range 0.97 to 1.04 for system sizes L = 10 to 60. 
There is no systematic variation with L and this indicates 
the absence of any finite size dependence for the above 
exponent at pc. Thus within our numerical accuracy we 
find that S{pc) — ^-0- It is thus clear from the Fig. |ll| 
that 6{pc) is the minimum of the p-dependent exponent 
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5{p). Note that the nonhnearity exponent a for the I-V 
curves would be just a = S + 1 = 2.0 (in this case), where 
/ ^ (V — Vg)", and for p sufficiently far from pc, a would 
also show the same concentration dependence. In com- 
parison, we observe that in an experiment on 2D arrays 
of normal metal islands connected by small tunnel junc- 
tions, Rimberg et al. found the nonlinearity exponent 
a to be 1.80 ± 0.16. It may also be noted that Roux 
and Herrmann | |2^ also obtained from their simulation 
of the I-V nonlinearity exponent a = 2 in their model 
2D network with the resistors without any dilution but 
with random thresholds. Clearly thus, our bond-diluted 
RRTN model is richer than the random threshold model 
at least as far as the nonlinearity exponent is concerned. 
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FIG. 11. The variation of nonlinearity exponent (5) for the 
G-V curves with p for a system of size 40 x 40. 

Now we discuss the crossover exponent which is an al- 
ternative way of accounting for the strength of nonlin- 
earity. The crossover exponent (cc) is defined from the 
power-law relationship Ic Gq where Ic is the crossover 
current at which the conductance of the system has in- 
creased from a nonzero initial value Go by a small but 
arbitrarily fixed fraction e. So x is defined only above Pc- 
The value of the above exponent was calculated by Gefen 
et. al., Q from the experimental nonlinear response data 
of discontinuous thick gold films near and above the per- 
colation threshold (in 3D). Their experimental measure- 
ment gives X = 1.47±0.10, whereas they argue through a 
model resistor network jl^] that the value of x should be 
3/2 (in 3D). The argument is based on the assumption of 
a power law dependence of conductance (G) which inter- 
polates between its initial value Gq (at I^ = V^) and the 
saturation value G/ {V = Vg). Note that the crossover 
exponent for the carbon-wax experiment in 3D is also 
close to this value. Gefen et al. |j] found that near the 
threshold S = t/v. Thus one expects that close to Pc, 
5 = 1.3/1.33 = 0.97 and the nonlinearity exponent for 
I-V characteristic is thus a = 5 + 1 = 1.97 in 2D. It will 
be noted that close to pc, the nonlinearity exponent for 
our model in 2D is very close to it. Further, since the 
exponent 5 lies between about 3 and 1 for p between 0.3 



and 0.9, the crossover exponent in our model can vary 
between 1.3 and 2 in the same dilution range in 2D. 

Apart from calculating the nonlinearity exponent as 
discussed previously, one may also measure or calcu- 
late the difference between the conductance (G) in the 
two asymptotic linear regimes. The two limiting conduc- 
tances Go and G / are already defined in two asymptotic 
linear regimes namely, V ^ Q and V oo respectively. 
Whereas the nonlinearity exponent [5 or a) acts as a 
measure of the initial nonlinearity near the threshold Vg , 
the difference in conductance Gd ~ G f ~ Gq serves as a 
measure of the overall nonlinearity. We have the follow- 
ing criteria which may be easily checked. 

• In the range p < pcf. Go = and G/ = 0, so 

Gd = 0, 

• In the range Pct < P < Pc- Go ~ and G/ 7^ 0, so 

Gd = Gf, 

• In the range Pc < p < 1: Gq > and Gd = 
Gf - Go. 
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FIG. 12. The behavior of Go against p. The system sizes 
(L) are indicated in the figure. 
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FIG. 13. The behavior of G/ against p. The system sizes 
(L) are indicated in the figure. 
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We have shown the variation of Go and G / against p in 
Fig. ^ and in Fig. |l^ respectively for a square network of 
size L — 12, 20, 28, 40 and 60. The averages are done over 
100 to 1000 configurations in each case. The interesting 
point to note is that as a function of p, G/ tends to be fiat 
(p-independent) as p \ but Go seems to be increasing 
sharply. Further, one may note that while Go behaves in 
a power-law fashion around pc = 0.5, G/ does the same 
around Pet — 0.181. Next, the difference in conductance 
Gd is plotted against p in Fig. |lj. The peak in this 
curve appears at around p ~ pc- This indicates that 
the overall nonlinearity is maximum near the geometrical 
percolation threshold. There does not seem to be any 
strong finite-size effects in the Gd-p graph except close 

to Pet- 




FIG. 14. The behavior of Gd against p for different system 
sizes L as indicated in the figure. 

The next important question is how the difference in 
conductance Gd is related to the initial conductance Go. 
The variation of Gd against Gq is shown in the Fig. |l^ 
for the whole range of p (0 to 1). As expected, the peak 
appears at around p = pct , and a strong finite size effect 
may be noted in the Gd-Go graph. One may also observe 
that in the thermodynamic limit, Gd = for p < pct. 
Now we would like to examine the above relationship, 
i.e., Gd against Go in the interval Pct < p < 1. We do 
fit the data with the function Gd ^ A — BGq{L), where 
A and B are constants and the exponent y{L) > 0. The 
data were fitted for different system sizes from i = 8 to 
80. In the Fig. |l^ we have shown y{L) against L to show 
the finite size behavior and find that y(oo) = 1.1 which 
is close to 1. Thus we may conclude that Gd is almost 
linearly dependent on Go in the limit of L ^ oo which 
means that G/ is also linearly dependent on Go- This 
in turn supports the idea of identical scaling relationship 
for two saturation conductances Go and G/ around the 
respective thresholds (pc and Pct). This is consistent with 
the fact that the system is in the same universality class 
in both the limits {i.e., Gf ~ (p — PdY). 




-O.l O.O O.l 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 



FIG. 15. The Gd against Go for different L as indicated in 
the figure. 
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FIG. 16. The finite size behavior of y{L), for different L 
data fitted with y{L) = 2/(00) + B / U where y{oo) = 1.1, B 
= 1.2, r = 0.5 



C. Comparison with some EM A results 

For the sake of comparison we plot the Gd against p 
and Gd against Gq curves with the corresponding results 
obtained by EMA for our model (in 2D). In Fig. |l^ the 
Gd against p is plotted for both the simulation result and 
the EMA result. Both are shown for our model network 
in 2D, where the simulation result plotted is for L = 60. 
In Fig. |l8|we plot Gd against Go again for the simulation 
and the EMA result. It is apparent from the above two 
figures that these results are close. 

The effective conductance of the network Gg = Gq 
when all the tunneling bonds behave as perfect insulators 
(the usual limit of binary mixture of conductors and insu- 
lators, Pc = 1/2) and the effective conductance Ge = G/ 
when all the tunneling bonds (or resistors) behave as the 
ohmic c ondu ctors (Eq. ( 3.1C )). Clearly, if one puts Pt = 
in Eq. ( [3.10 ), one finds that 
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Go — 



dPg-l 

id-l)- 



(4.11) 



The behavior of Gd against p is shown in Fig. 17 for 2D. 
Gd has a maximum at pc (= 1/2 in 2D) as it should 
be because the system with the o- bonds only is most 
tenuous at the geometrical percolation threshold and the 
measure of nonlinearity clearly should be maximum there 
(with the largest concentration of the t-bonds). The plot 
of Gd against Gq is also shown in Fig. In the same plot 
we have also shown the simulation result for comaprison 
and one may note that the quantity Gd is maximum at 
p = Pc- The EMA result is quite close to the simulation 
result and the match becomes better as one goes farther 
and farther away from pd- 




FIG. 17. The variation of Gd against p, obtained by the 
simulation and the EMA, are plotted for comparison. 
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FIG. 18. The variation of Gd against Go, obtained by the 
simulation and the EMA, are plotted for comparison. 



V. SUMMARY AND COMMENTS 

In this paper, we have been mainly concerned with 
the nonlinear response chracteristics, the relevant scaling 



and related exponents in composites or granular metallic 
systems where transport due to charge tunneling plays 
an important role. One may note that our correlated 
percolation model and the corresponding RRTN network 
can capture the essential features related to nonlinear- 
ity quite well. Before ending we would like to make 
some comments and discussions for the sake of complete- 
ness. The existence of a well-defined voltage threshold 
(Vg) close to which a power-law regime appears indi- 
cates that this threshold (breakdown) voltage is a critical 
point. Similar threshold behavior has been observed in 
pinned charge-density-wave (CDW) conductors, super- 
conductors with pinned vortices (type-II), etc. Indeed 
it has been argued long ago that wherever there exists 
some sort of threshold of force for motion to occur, the 
threshold actually corresponds to a dynamic critical point 
p5| for the driven dynamical system. Disorder in such 
systems is known to give rise to 'pinning' or inhibition 
to transport upto a critical value of the force. Clearly in 
our percolation model, the threshold Vg acts as a dynam- 
ical critical point for the systems with volume fractions 
in the range Pct < p < Pc- For such systems the field 
corresponding to the threshold voltage Vg is called the 
dielectric breakdown field and is pretty well-studied in 
the DRRN type models. We have focussed on the 
threshold as a dynamical critical point in our discrete 
model, and calculated the breakdown exponent for our 
RRTN model elsewhere [!§. 

As we pointed out above (see also Ref. [Q) that the 
mechanism of nonlinearity is essentially the same both 
below and above the system threshold. Below the sys- 
tem thresold, there is no system spanning cluster. So the 
transport is identified as intercluster tunneling or hop- 
ping across dangling bonds or gaps. Above the thresh- 
old there are both intercluster and intracluster tunneling. 
But intracluster tunneling mechanism certainly domi- 
nates. The nearest neighbor gaps are everywhere; both 
inside the smaller isolated clusters as well as in the system 
spanning cluster. So the tunneling mechanism is opera- 
tive both below and above Pc (in the interval Pct < P < 1) 
giving rise to nonlinear regime in the response. 

From our analysis of current-voltage [I-V) data we 
understand that one does not do justice by fitting only 
the I-V curve and finding out the nonlinearity exponent 
therefrom since that fitting is not robust. One may eas- 
ily get tempted to fit the nonlinear regime of a I-V curve 
in general through an n-th degree polynomial function. 
A reasonable choice Q would be to fit with the law: 
/ = GiV -I- G2V^, assuming that the leading nonlinear 
term is cubic (ignoring other higher order terms). Even 
ignoring the fact that the nonlinearity exponent away 
from Pc are significantly different from integer values, we 
note that this type of analysis in this form may lead to 
confusion. The determination of the exponents and the 
coefficients may become arbitrary. Arbitrary selection 
of the range of I-V data for this purpose and the fit- 
ting of that would account for this arbitrariness. One 
may not know upto what voltage scale (in the nonlin- 
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ear regime) one should fit the data and hence one may 
get a set of different answers for the values of the expo- 
nents. A typical I-V curve, for the kind of systems we 
address, can in general be fitted by a simple power-law: 
J (Y — Vg)"-, at least around the threshold voltage 
Vg. But this also may not provide unambiguous answer 
most of the time. One may note that the exponent for 
such power-law does change continuously as the applied 
voltage V is increased from Vg and the I-V curve may 
ultimately approach another linear regime. In fact, the 
selection of the range of data and its fitting had been at 
the root of many confusing results as found in the lit- 
erature. Instead our prescription of fitting of the G-V 
data for the entire nonlinear regime along with the sat- 
uration, as presented, is fo und to be most satisfactory. 
^From such a fitting (Eq. (4.7)) of the G-V data, one 
can obtain the desired power-law and find the exponent. 
The G-V curves for the type of composites we focuss our 
attention on, may all be generically fitted by a function 
like G{V) ^Gn + Gdf{AV), where f{AV) is a function 
that behaves as AV^'' for small AV, where AV ^V -Vg. 
f{AV) as AV^ ^ and f{AV) ^ 1 as Ay ^ oo (or 
an appropriately large value for a given finite system). 

Next we point out a few observations on some recent 
experiments in the literature which are found to be not 
in full agreement with our simulation results. In a re- 
cent experiment on the carbon-wax mixture [ p7| it is 
claimed that the conductivity exponent {t) to be dif- 
ferent in the two extreme limits namely, AV^ — > and 
AV oo. More explicitly and in the perspective of the 
foregoing discussion, Go ^ Ap* and Gf ^ Ap* , where 
Ap = (p — pc). Note that the Ap for the scaling of G/ 
should have been {p — pct) as we have discussed before. 
In any case, it is claimed that t' ^ t. More specifically, 
t' = ct, where c 0.76 (in 3D) as reported This 
in turn is related to the question of universality class, 
claiming that the system goes over to different univer- 
sality class in the saturation limit [V —>■ oo ) ( has been 
termed as altered percolatin g stat e in Ref. |2^). How- 
ever, ws we discuss in section [VB (we have also claimed 
earlier |Q ) that the system seems to be in the same uni- 
versality class in both the limits whereby we claim that 
t = t' within the numerical accuracy. Our calculation of 
course is in 2D. But the essential physics should remain 
the same as we go over to 3D. 

We have seen that the nonlinearity exponent, a (= 
(5-1-1) varies significantly with the volume fraction (p) of 
the conducting component, the minimum value of the ex- 
ponent being around 2.0 at p = pc- Chen and Johnson, in 
their experiment on Ag-KCl (t) , found the above nonlin- 
earity exponent to vary with p. In particular they found 
that for a sample very close to percolation threshold, the 
nonlinearity exponent is as large as 20. In another early 
and a very prominent experiment, nonlinear response in 
ZnO varistors has been addressed by Mahan et al. [ p9[ . 
The observed I-V characteristics in the ZnO varistors are 
often empirically described by the power-law relation: 



/ = kV". 



(5.1) 



where the parameter a (> 1) is the measure of nonlin- 
earity. A theory had been developed by Mahan et al. to 
predict the coefficient of nonlinearity (a) as high as 50 or 
even 100. Such a high value of a of course indicates an ex- 
ponential relationship rather than a power-law. Canessa 
and Nguyen |Q had attempted a computer simulation 
study of the varistors considering a binary mixture model 
to understand the very high value of a. 

Intrigued by these results we ran some test for a square 
lattice with (L=20 and 40) p=0.2, very close to our pct. 
We obtained in our preliminary analysis by fitting G with 
voltage V that a > 20 for L = 20 and a ^ 1 (diverging 
for ever) for L = 40. This was very suspicious and by 
zooming into the fitting close to the threshold (Vg), we 
found the fitting to be extremely bad. The other func- 
tions listed before (including the double-exponential) did 
not give any better fitting either. But this is the region 
which gives the initial power-law exponent. To get the 
fit better in this region, we fitted the logarithmic con- 
ductance (In AG), keeping our fitting function the same 
Eq. (4/7) as before. Then we found the exponent 5 to be 
of the order of 3 (the Fig. [ill was done by taking care 
of this fitting problem). Thus power-law description still 
seems to hold but there seems to be a different problem 
close to Pct- Since the fitting was still not as good as at 
higher p > 0.4, we wanted to check if the averaging pro- 
cess itself is at flaw at very low volume fractions. Hence 
we look at the distribution of the current (/) for different 
realizations of the sample of size L = 20 and p ~ 0.2 for 
two different voltages. In Fig.|l^, the histogram is shown 
for y = 5, and in Fig.^ the same is shown for V — 10. 
The distribution for y=10 is reasonably well behaved, 
but the one for V=5 has an isolated delta-function-like 
peak at zero conductance apart from two other broader 
peaks at higher conductances. We did also look at the rel- 
ative value of the variance (relative to the mean squared) 
defined as 



RV = 



(<P>-<I >^) 
< / >2 



(5.2) 



If this value is less than 1, averaging is alright, but if it 
is much larger than 1, the self- averaging property of the 
distribution does not work. We have tested that RV is 
about 3.3 for the set of data for y = 5 and it is about 0.3 
for y=10. Thus the low p samples are non-self-averaging 
at low voltages but tend to be more self-averaging at 
higher voltages. This property was also reported for 
quantum systems by Lenstra and Smokers |^l| , and thus 
our semi-classical model of percolation indeed seems to 
capture some non-classical (quantum) behavior at low p 
as expected. 

The connection between the above exponent 5 and the 
crossover exponent {x) is the following. The crossover 
voltage Vc is related to the initial nonzero conductance 
Go by Vc ~ Go"^ But from Eq. ( |4.10| ) we know that 
the excess conductance AG(y) = Ay" for small Ay. 
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Choosing AG = eGo for an arbitrary but small e and if 
one is close to percolation threshold, one has Vc ~ Gq^^ 
which implies that 5 = l/{x — 1). Thus, as the value of 
nonlinearity exponent is dependent on p, the crossover 
exponent (x) should also be dependent on p. The values 
of the crossover exponent x so far reported in the litera- 
ture are for the experimental data for the samples close 
to but above the threshold. But we may conclude that 
this exponent should also show change with p if p is away 
frompc- In fact, such an example may be mentioned here. 
The crossover exponent is found to be widely different for 
a set of samples as measured recently ||2^. The reason 
for this is yet to be understood. But if the samples taken 
for measurements have widely different volume fractions 
(p) of conducting components, then nonlinearity expo- 
nent and thus the crossover exponent for them could be 
widely different. 




Now it may be note d that if 5 > 1 or in other words 
/i7 > 1 in the Eq. (4/7), the first derivative of G with re- 
spect to V, i.e., dG/dV would attain a maximum value 
at the inflexion point (at some voltage in the nonlinear 
regime). Wc fitted a set of G-V and the corresponding 
dG/dV-V data (see Ref. ^^), taken from the experi- 
mental observations on a carbon-wax composite system 
to show the appearance of a peak in the dG/dV-V 
curve. To indicate how good or bad our fitting equation 
is, we fit the above data set for G-V characteristics by the 
Eq. ( [4.7] ). The fitted line is seen to match the experimen- 
tal data very well with the above function. As a further 
test of the fitting, we took the parameters of the G-V^-fit, 
and used them (the functional form) to obtain dG/dV as 
a function of V. The agreement with the experimental 
values for dG/dV should be considered rather good given 
that G and dG/dV are independent measurements and 
that higher the harmonic more error-prone is its value. 
Incidentally, the nonlinearity exponent 5 for this 3D ex- 
perimental data fitted by our method comes out to be 
1.74, where the crossover exponent measurement on the 
same sample would give an 5- value of about 2. 
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FIG. 19. The distribution of current (/) in the 2D square 
network of size L = 20 for p = 0.2 and an applied voltage V 
— 5.0. This distribution is non-self-averaging. 6x10* config- 
urations were used for this purpose. 




FIG. 20. The distribution of current (/) in the 2D square 
network of size L — 20 for p = 0.2 and an applied voltage V — 
10.0. This distribution is self-averaging. 2 x 10'* configurations 
were used for this purpose. 
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